/*
 *  linkage.c
 *  Sandbox
 *
 *  Created by Joshua Lewis on 10/28/10.
 *  A parallel single linkage algorithm using OpenMP.
 *
 */

#include "linkage.h"

int compareFloat(const void *a, const void *b) {
	float a_float = *(float *)a;
	float b_float = *(float *)b;
	
	if (a_float < b_float) return -1;
	if (a_float > b_float) return 1;
	return 0;
}

int compareDendrite(const void *a, const void *b) {
	float *a_float = &(*(dendrite *)a).distance;
	float *b_float = &(*(dendrite *)b).distance;
	
	return compareFloat(a_float, b_float);
}

void nodist_dendrogram(int N, int remainSize, dendrite *sorted, dendrite *result) {
	int i, j;
	int totalSize = 0;
	int baseSize = 0;
	int pruneSize = 0;
	int bookmark = 0;
	int stop = 0;
	int threadNum = omp_get_max_threads(); // Number of threads OpenMP will spawn
	int branchSize = 500;
	
	dendrite * branch = (dendrite *)malloc(threadNum * (branchSize + N - 1) * sizeof(dendrite));
	dendrite * temp = (dendrite *)malloc(threadNum * branchSize * sizeof(dendrite));
	int * count = (int *)calloc(threadNum, sizeof(int));
	
#pragma omp parallel private(i, j)
	{
		int th_id = omp_get_thread_num();
		int branch_index = th_id * (branchSize + N - 1);
		int temp_index = th_id * branchSize;
		int sortedStart = 0;
		
		while (stop != 1 && remainSize - bookmark + pruneSize >= threadNum * branchSize) {
			
			for (i = 0; i < baseSize; i++) //Copy in base
				branch[branch_index + i] = sorted[i];
			
			for (i = 0; i < branchSize; i++) // Copy in branch
			{
				if (i + (th_id * branchSize) < pruneSize)
					sortedStart = baseSize;
				else
					sortedStart = bookmark - pruneSize; // Since we want to start at bookmark and i + (th_id * branchSize) >= pruneSize
				
				branch[branch_index + baseSize + i] = sorted[sortedStart + (th_id * branchSize) + i];
				temp[temp_index + i] = sorted[sortedStart + (th_id * branchSize) + i];
			}
			
			count[th_id] = prune(&branch[branch_index], N, branchSize + baseSize) - baseSize;
			
			
#pragma omp barrier  
			
			int offset = 0; // Offset in memory based on higher threads' counts
			for(i = 0; i < th_id; i++)
				offset += count[i];
			
			j=0;
			for (i = 0; i < branchSize; i++) // Copy real joins that aren't from the base
				if(branch[branch_index + baseSize + i].i != -1) {
					sorted[baseSize + offset + j] = temp[temp_index + i];
					j++;
				}
			
#pragma omp barrier
			
			if(th_id == 0) {
				bookmark += threadNum * branchSize - pruneSize; // Move bookmark forward by the number of elements we read off
				totalSize = 0;
				for(i = 0; i < threadNum; i++)
					totalSize += count[i];
				baseSize += count[0];
				pruneSize = totalSize - count[0];
				if(baseSize == N - 1) // If we're done, stop
					stop = 1;
			}
			
#pragma omp barrier
		}
	}
	
	
	if(baseSize < N - 1)
		for(i = 0; i < remainSize - bookmark; i++) // Copy over remainder
			sorted[baseSize + pruneSize + i] = sorted[bookmark + i];
	else
		pruneSize = remainSize = bookmark = 0; // Just prune the base to get group numbers
	
	
	prune(sorted, N, baseSize + pruneSize + remainSize - bookmark);
	
	j = 0;
	for (i = 0; i < baseSize + pruneSize + remainSize - bookmark; i++)
		if (sorted[i].i > -1)
			result[j++] = sorted[i];
	
	free(sorted);
	free(branch);
	
	return;
}

void comp(int N, float *distance, dendrite *result) {
	int i, j, k;
	int * nearest = (int *)malloc((N - 1) * sizeof(int));
	int * bookkeep = (int *)malloc(N * sizeof(int));
	int * ncopy = (int *)malloc((N - 1) * sizeof(int));
	float * minima = (float *)malloc((N - 1) * sizeof(float));
	float min = FLT_MAX;
	float ndist;
	int mindex, groupA, groupB;
	
	//step 1: make two arrays of size n-1: row minima and nearest neighbors
	for(i = 0; i < N - 1; i++) {
		min = FLT_MAX;
		for(j = i + 1; j < N; j++) 
			if (distance[j*(j-1)/2+i] < min) {
				min = distance[j*(j-1)/2 +i];
				mindex = j;
			}
		minima[i] = min;
		nearest[i] = mindex;
		ncopy[i] = mindex;
	}
	
	for(i = 0; i < N; i++)
		bookkeep[i] = i;
	
	//	for(i = 0; i < (N-1); i++)
	//		printf("%i \t", i);
	//	printf("\t\t\t");
	//	for(i = 0; i < N; i++)
	//		printf("%i \t", i);
	//	printf("\t\t\t");
	//	for(i = 0; i < (N-1); i++)
	//		printf("%i \t", i);
	//	printf("\n\n");
	
	// steps 2-5
	for (j = 0; j < (N-1); j++) {
		//for (i = 0; i < N - 1; i++)
		//printf("%i %f \n", nearest[i], minima[i]);
		min = FLT_MAX;
		//find minimum and index of minimum in minima
		for (i = 0; i < (N-1); i++) {
			if (minima[i] < min) {
				min = minima[i];
				mindex = i;
			}
		}
		
		result[j].i = bookkeep[mindex];
		result[j].j = nearest[mindex];
		result[j].distance = min;
		
		if(nearest[mindex] > N - 1) // check for group origin and switch if necessary
			for(i = 0; i < N - 1; i++) {
				if(bookkeep[i] == nearest[mindex] && i < mindex) {
					result[j].i = result[j].j;
					result[j].j = bookkeep[mindex];
					ncopy[i] = mindex;
					mindex = i;
				}
				if(bookkeep[i] == nearest[mindex] && i < ncopy[mindex])
					ncopy[mindex] = i;
			}
		
		//		for(i = 0; i < (N-1); i++)
		//			printf("%i \t", nearest[i]);
		//		printf("\t\t\t");
		//		for(i = 0; i < N; i++)
		//			printf("%i \t", bookkeep[i]);
		//		printf("\t\t\t");
		//		for(i = 0; i < (N-1); i++)
		//			printf("%i \t", minima[i] == FLT_MAX);
		//		printf("\n");
		
		groupA = bookkeep[mindex];
		groupB = bookkeep[ncopy[mindex]];
		
		//absorb cluster and update nearest neighbors
		for (i = 0; i < (N-1); i++) {
			//	printf("%i %i %i \n", nearest[i], ncopy[mindex], bookkeep[mindex]);
			if (nearest[i] == groupA || nearest[i] == groupB && i != mindex && i != ncopy[mindex])
				nearest[i] = N + j;				
		}
		
		for(i = 0; i < N; i++) {
			if (bookkeep[i] == groupA || bookkeep[i] == groupB)
				bookkeep[i] = N + j;
			if (i > ncopy[mindex]) {
				ndist = distance[i*(i-1)/2 + ncopy[mindex]];
				if (ndist < distance[i*(i-1)/2 + mindex])
					distance[i*(i-1)/2 + mindex] = ndist;
			}
			if(i > mindex && i < ncopy[mindex]) {
				ndist = fmin(distance[i*(i-1)/2 + mindex], distance[ncopy[mindex] * (ncopy[mindex] - 1) / 2 + i]);
				distance[i*(i-1)/2 + mindex] = ndist;
				//distance[ncopy[mindex] * (ncopy[mindex] - 1) / 2 + i] = ndist;
			}
			if(i < mindex && i < ncopy[mindex]) {
				ndist = fmin(distance[mindex * (mindex - 1) / 2 + i], distance[ncopy[mindex] * (ncopy[mindex] - 1) / 2 + i]);
				distance[mindex * (mindex - 1) / 2 + i] = ndist;
				//distance[ncopy[mindex] * (ncopy[mindex] - 1) / 2 + i] = ndist;
			}
		}
		
		//replace minima[NN[i]] with FLT_MAX
		minima[ncopy[mindex]] = FLT_MAX;
		distance[ncopy[mindex]*(ncopy[mindex]-1)/2 + mindex] = FLT_MAX;
		
		//		for (i = 0; i < N-1; i++)
		//			if (bookkeep[i] == N + j && i != mindex)
		//				minima[i] = FLT_MAX;
		
		//find new minimum of row i
		i = mindex;
		min = FLT_MAX;
		for (k = i + 1; k < N; k++)
			if (distance[k*(k-1)/2 + i] < min && minima[k] != FLT_MAX) {
				min = distance[k*(k-1)/2 +i];
				mindex = k;
				//printf("%i %i %i \t",i,  k, ncopy[i]);
			}
		//printf("\n");
		minima[i] = min;
		nearest[i] = bookkeep[mindex];
		ncopy[i] = mindex;
		
	}		
	
	free(minima);
	free(ncopy);
	free(nearest);
	free(bookkeep);
	return;
}


void dendrogram(int N, float *distance, dendrite *result) {
	int i, j;
	int totalSize = 0;
	int baseSize = 0;
	int pruneSize = 0;
	int bookmark = 0;
	int stop = 0;
	int remainSize = (N * (N - 1) / 2);
	int threadNum = omp_get_max_threads(); // Number of threads OpenMP will spawn
	int branchSize = 500;
	
	dendrite * sorted = (dendrite *)malloc((N * (N - 1) / 2) * sizeof(dendrite));
	dendrite * branch = (dendrite *)malloc(threadNum * (branchSize + N - 1) * sizeof(dendrite));
	dendrite * temp = (dendrite *)malloc(threadNum * branchSize * sizeof(dendrite));
	int * count = (int *)calloc(threadNum, sizeof(int));
	
	for(i = 0; i < N - 1; i++)
		for(j = i + 1; j < N; j++) {
			sorted[j * (j - 1) / 2 + i].distance = distance[j * (j - 1) / 2 + i];
			sorted[j * (j - 1) / 2 + i].i = i;
			sorted[j * (j - 1) / 2 + i].j = j;
		}
	
	psort(sorted, N * (N - 1) / 2, sizeof(dendrite), compareDendrite);
	
#pragma omp parallel private(i, j)
	{
		int th_id = omp_get_thread_num();
		int branch_index = th_id * (branchSize + N - 1);
		int temp_index = th_id * branchSize;
		int sortedStart = 0;
		
		while (stop != 1 && remainSize - bookmark + pruneSize >= threadNum * branchSize) {
			
			for (i = 0; i < baseSize; i++) //Copy in base
				branch[branch_index + i] = sorted[i];
			
			for (i = 0; i < branchSize; i++) // Copy in branch
			{
				if (i + (th_id * branchSize) < pruneSize)
					sortedStart = baseSize;
				else
					sortedStart = bookmark - pruneSize; // Since we want to start at bookmark and i + (th_id * branchSize) >= pruneSize
				
				branch[branch_index + baseSize + i] = sorted[sortedStart + (th_id * branchSize) + i];
				temp[temp_index + i] = sorted[sortedStart + (th_id * branchSize) + i];
			}
			
			count[th_id] = prune(&branch[branch_index], N, branchSize + baseSize) - baseSize;
			
			
#pragma omp barrier  
			
			int offset = 0; // Offset in memory based on higher threads' counts
			for(i = 0; i < th_id; i++)
				offset += count[i];
			
			j=0;
			for (i = 0; i < branchSize; i++) // Copy real joins that aren't from the base
				if(branch[branch_index + baseSize + i].i != -1) {
					sorted[baseSize + offset + j] = temp[temp_index + i];
					j++;
				}
			
#pragma omp barrier
			
			if(th_id == 0) {
				bookmark += threadNum * branchSize - pruneSize; // Move bookmark forward by the number of elements we read off
				totalSize = 0;
				for(i = 0; i < threadNum; i++)
					totalSize += count[i];
				baseSize += count[0];
				pruneSize = totalSize - count[0];
				if(baseSize == N - 1) // If we're done, stop
					stop = 1;
			}
			
#pragma omp barrier
		}
	}
	
	
	if(baseSize < N - 1)
		for(i = 0; i < remainSize - bookmark; i++) // Copy over remainder
			sorted[baseSize + pruneSize + i] = sorted[bookmark + i];
	else
		pruneSize = remainSize = bookmark = 0; // Just prune the base to get group numbers
	
	
	prune(sorted, N, baseSize + pruneSize + remainSize - bookmark);
	
	j = 0;
	for (i = 0; i < baseSize + pruneSize + remainSize - bookmark; i++)
		if (sorted[i].i > -1)
			result[j++] = sorted[i];
	
	free(sorted);
	free(branch);
	
	return;
}

int prune (dendrite *branch, int N, int size) {
	int i, j;
	int groupCount = 0;
	
	for(i = 0; i < size; i++)
	{
		if (branch[i].i > -1)
			if (branch[i].i == branch[i].j)
				branch[i].i = -1;
			else {
				for(j = i + 1; j < size; j++) {
					if(branch[j].i > -1) {
						if ((branch[i].i == branch[j].i) || (branch[i].j == branch[j].i))
							branch[j].i = N + groupCount;
						if ((branch[i].i == branch[j].j) || (branch[i].j == branch[j].j))
							branch[j].j = N + groupCount;
					}
				}
				groupCount++;
			}
	}
	
	return groupCount;
}

void assignLaunch(dendrite *dendrogram, int k, int N, int *result) {
	int curK = 0;
	for(int i = 0; i < k - 1; i++)
	{
		if(dendrogram[N - i - 2].i < N)
			result[dendrogram[N - i - 2].i] = curK++;
		else if(dendrogram[N - i - 2].i <= 2 * (N - 1) - (k - 1))
			assign(dendrogram, i - 2, curK++, N, result);
		
		if(dendrogram[N - i - 2].j < N)
			result[dendrogram[N - i - 2].j] = curK++;
		else if(dendrogram[N - i - 2].j <= 2 * (N - 1) - (k - 1))
			assign(dendrogram, i - 2, curK++, N, result);    
	}
}

void assign(dendrite *dendrogram, int line, int k, int N, int *result) {
	if(dendrogram[line].i < N)
		result[dendrogram[line].i] = k;
	else
		assign(dendrogram, dendrogram[line].i - N, k, N, result);
	
    if(dendrogram[line].j < N)
		result[dendrogram[line].j] = k;
    else
		assign(dendrogram, dendrogram[line].j - N, k, N, result);
}

void denSlow(int N, float *distance, dendrite *result) {
	int i, j;
	
	int *tracker = (int *)malloc(sizeof(int) * N);
	float *index = (float *)malloc(sizeof(float) * N * (N - 1));
	int *dend = (int *)malloc(sizeof(int) * 2 * (N - 1));
	
	
	float min;
	int hold;
	int a, b, akIndex, bkIndex;
	
	for(i = 0; i < N; i++)
	{
		tracker[i] = i;
		for(j = i + 1; j < N; j++)
		{
			index[j * (j - 1) + 2 * i] = i;
			index[j * (j - 1) + 2 * i + 1] = j;
		}
	}
	
	
	for(j = 0; j < N - 1; j++) {
		min = MAXFLOAT;
		for(i = 0; i < N * (N - 1) / 2; i++) { // Find minimum distance
			if (distance[i] < min) {
				min = distance[i];
				hold = i;
			}
		}
		
		a = index[2 * hold];
		b = index[2 * hold + 1];
		
		dend[2 * j] = tracker[a];
		dend[2 * j + 1] = tracker[b];
		//dend[3 * j + 2] = min;
		
		result[j].i = tracker[a];
		result[j].j = tracker[b];
		result[j].distance = min;
		
		distance[hold] = MAXFLOAT;
		
		for (i = 0; i < N; i++) { // Update distance matrix
			
			if (a > i)
				akIndex = a * (a - 1) / 2 + i;
			else
				akIndex = i * (i - 1) / 2 + a;
			
			if (b > i)
				bkIndex = b * (b - 1) / 2 + i;
			else
				bkIndex = i * (i - 1) / 2 + b;
			
			
			if (i == a)
				tracker[i] = N + j;
			else if (i == b)
				tracker[i] = -1;
			else if ((tracker[i] != -1) && (i != a) && (i != b)) {
				distance[akIndex] = fmin(distance[akIndex], distance[bkIndex]);
				distance[bkIndex] = MAXFLOAT;
			}
		}
	}
  	
	free(tracker);
	free(index);
	free(dend);
	
	return;	
	
}

void denSlow2(int N, float *distance, dendrite *result) {
	int i, j;
	int denIndex = 0;
	int remainSize = (N * (N - 1) / 2);
	int threadNum = omp_get_max_threads(); // Number of threads OpenMP will spawn
	int branchSize = remainSize / threadNum;
	
	dendrite * sorted = (dendrite *)malloc((N * (N - 1) / 2) * sizeof(dendrite));
	dendrite * branch = (dendrite *)malloc(threadNum * branchSize * sizeof(dendrite));
	dendrite * temp = (dendrite *)malloc(threadNum * branchSize * sizeof(dendrite));
	int * count = (int *)malloc(threadNum * sizeof(int));
	
	for(i = 0; i < N - 1; i++)
		for(j = i + 1; j < N; j++) {
			sorted[j * (j - 1) / 2 + i].distance = distance[j * (j - 1) / 2 + i];
			sorted[j * (j - 1) / 2 + i].i = i;
			sorted[j * (j - 1) / 2 + i].j = j;
		}
	
	psort(sorted, N * (N - 1) / 2, sizeof(dendrite), compareDendrite);
	
#pragma omp parallel private(i)
	{
		int th_id = omp_get_thread_num();
		int th_index = th_id * branchSize;
		
		for (i = 0; i < branchSize; i++)
		{
			branch[th_index + i] = sorted[(branchSize * th_id) + i];
			temp[th_index + i] = sorted[(branchSize * th_id) + i];
			//printf("%d %d %f\n", branch[th_id * branchSize + i].i, branch[th_id * branchSize + i].j, branch[th_id * branchSize + i].distance);
		}
		
		count[th_id] = prune(&branch[th_index], N, branchSize);
		
#pragma omp barrier  
		
		int offset = 0;
		for(i = 0; i < th_id; i++)
			offset += count[i];
		
		for (i = 0; i < branchSize; i++)
			if(branch[th_index + i].i != -1)
				sorted[offset + i] = temp[th_index + i];
	}
	
	int sum = 0;
	for(i = 0; i < threadNum; i++)
		sum += count[i];
	
	int leftover = remainSize % branchSize;
	for(i = 0; i < leftover; i++)
		sorted[sum + i] = sorted[branchSize * threadNum + i];
	
	prune(sorted, N, sum + leftover);
	
	for (i = 0; i < sum + leftover; i++)
		if (sorted[i].i > -1) {
			result[denIndex] = sorted[i];
			denIndex++;
		}
	
	//for(i = 0; i < N - 1; i++)
	//  printf("%d %d %f\n", result[i].i, result[i].j, result[i].distance);
	
	
	free(sorted);
	free(branch);
	
	return;
}
